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DNS  and  LES  of  a shear-free  mixing  layer 

By  B.  Knaepen,  0.  Debliquyt  and  D.  Carati  f 


1.  Introduction 

The  shear-free  mixing  layer  represents  one  of  the  simplest  inhomogeneous  flows.  It 
consists  of  two  ‘distinct’  homogeneous  regions  of  different  turbulent  kinetic  energy  in- 
teracting through  a layer  of  rapid  transition.  The  layer  is  said  to  be  shear-free  since  the 
two  homogeneous  regions  have  no  relative  velocity.  While  flows  encountered  in  nature 
or  industrial  applications  are  more  often  not  devoid  of  shear,  the  study  of  the  shear-free 
mixing  layer  is  nevertheless  useful  since  it  allows  the  mixing  properties  of  turbulence  to 
be  examined  in  a simplified  environment.  Indeed,  turbulent  transport  properties  in  shear 
flows  are  more  difficult  to  track  since  they  can  be  overwhelmed  by  production  sources 
originating  from  gradients  in  the  mean  velocity. 

The  shear-free  mixing  layer  has  already  received  attention  in  the  past,  both  from  the 
experimental  and  the  numerical  point  of  view.  The  first  experimental  study  of  the  shear- 
free  mixing  layer  is  due  to  Glibert  (1980).  The  flow  was  obtained  by  forcing  a stream 
through  a grid  with  two  different  mesh  spacings.  The  two  sides  of  the  grid  have  however 
equal  solidity  resulting  in  an  outgoing  shear-free  flow.  In  this  first  experimental  study,  the 
author  mainly  concentrated  his  attention  on  the  downstream  evolution  of  the  spreading- 
rate  parameter,  which  is  a measure  of  the  thickness  of  the  mixing-layer.  A more  extensive 
experimental  study  was  later  performed  by  Veeravalli  & Warhaft  (1987,  1989).  Owing  to 
a different  experimental  setup,  the  authors  achieved  a higher  ratio  between  the  energies 
characterizing  the  two  sides  of  the  flow  and  this  resulted  in  the  observation  of  large-scale 
intermittency  in  the  flow.  The  Veeravalli  & Warhaft  (1989)  data  will  be  used  here  as  the 
main  benchmark  for  the  present  study  since  it  contains  a detailed  documentation  of  the 
flow  characteristics  we  will  be  examining.  Prom  the  numerical  point  of  view,  the  shear-free 
mixing  layer  was  studied  in  Briggs  et  al.  (1996)  using  DNS.  In  that  article  the  authors 
also  used  the  data  of  Veeravalli  & Warhaft  (1989)  as  the  point  of  comparison.  Their 
simulations  were  done  using  a spectral  code  with  a resolution  of  1283  Fourier  modes. 
The  microscale  Reynolds  numbers  reached  in  the  low-  and  high-energy  homogeneous 
regions  were,  respectively,  11  and  22.5,  which  are  roughly  half  of  the  values  reported 
in  the  Veeravalli  & Warhaft  (1989)  experiment  (note  that  Briggs  et  al.  (1996)  use  a 
different  definition  of  the  microscale  Reynolds  number  than  the  one  in  Veeravalli  & 
Warhaft  (1989)  and  here).  However,  when  properly  nondimensionalized,  they  were  able 
to  reproduce  satisfactorily  the  turbulence  statistics  of  the  flow. 

The  purpose  of  this  work  is  twofold.  First,  given  the  computational  resources  available 
today,  it  is  possible  to  reach,  using  DNS,  higher  Reynolds  numbers  than  in  Briggs  et  al. 
(1996).  In  the  present  study,  the  microscale  Reynolds  numbers  reached  in  the  low-  and 
high-energy  homogeneous  regions  are,  respectively,  32  and  69.  The  results  reported  earlier 
can  thus  be  complemented  and  their  robustness  in  the  presence  of  increased  turbulence 
studied.  The  second  aim  of  this  work  is  to  perform  a detailed  and  documented  LES  of 
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Figure  1.  Graphical  representation  of  the  shear-free  mixing  layer 


the  shear-free  mixing  layer.  In  that  respect,  the  creation  of  a DNS  database  at  higher 
Reynolds  number  is  necessary  in  order  to  make  meaningful  LES  assessments.  Prom  the 
point  of  view  of  LES,  the  shear- free  mixing-layer  is  interesting  since  it  allows  one  to 
test  how  traditional  LES  models  perform  in  the  presence  of  an  inhomogeneity  without 
having  to  deal  with  difficult  numerical  issues.  Indeed,  as  argued  in  Briggs  et  al.  (1996), 
it  is  possible  to  use  a spectral  code  to  study  the  shear-free  mixing  layer  and  one  can  thus 
focus  on  the  accuracy  of  the  modelling  while  avoiding  contamination  of  the  results  by 
commutation  errors  etc. 

This  paper  is  organized  as  follows.  First  we  detail  the  initialization  procedure  used  in 
the  simulation.  Since  the  flow  is  not  statistically  stationary,  this  initialization  procedure 
has  a fairly  strong  influence  on  the  evolution.  Although  we  will  focus  here  on  the  shear- 
free  mixing  layer,  the  method  proposed  in  the  present  work  can  easily  be  used  for  other 
flows  with  one  inhomogeneous  direction.  The  next  section  of  the  article  is  devoted  to 
the  description  of  the  DNS.  All  the  relevant  parameters  are  listed  and  comparison  with 
the  Veeravalli  & Warhaft  (1989)  experiment  is  performed.  The  section  on  the  LES  of  the 
shear-free  mixing  layer  follows.  A detailed  comparison  between  the  filtered  DNS  data  and 
the  LES  predictions  is  presented.  It  is  shown  that  simple  eddy  viscosity  models  perform 
very  well  for  the  present  test  case,  most  probably  because  the  flow  seems  to  be  almost 
isotropic  in  the  small-scale  range  that  is  not  resolved  by  the  LES. 


2.  Initialization  of  the  flow 

Prom  the  numerical  point  of  view,  one  of  the  most  appealing  properties  of  the  shear- 
free  mixing  layer  is  the  possibility  of  simulating  this  flow  with  a purely  spectral  three 
dimensional  code.  Indeed,  periodicity  can  be  enforced  by  considering  a second  mixing- 
layer,  which  performs  the  ‘reverse’  transition  compared  to  the  first  one.  The  situation  is 
depicted  in  Fig.  1.  This  also  has  the  advantage  that  results  gathered  from  the  two  mixing 
layers  can  be  averaged  to  improve  the  statistics.  This  possibility  will  be  systematically 
exploited  in  the  presentation  of  numerical  results  in  the  following  sections. 

Given  a 3D  spectral  code,  the  non- trivial  part  of  the  simulation  is  then  to  build  a 
suitable  initial  condition  that  mimics  the  mixing  layer.  Indeed,  because  the  decaying 
mixing  layer  is  not  statistically  stationary,  it  is  not  acceptable  to  wait  until  the  initial 
condition  is  forgotten  by  the  flow  due  to  the  stochastic  nature  of  turbulence.  Indeed, 
the  simulation  will  remain  quite  strongly  influenced  by  the  initial  state  of  the  velocity 
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field.  To  proceed,  it  is  necessary  to  introduce  a few  definitions  and  notations.  The  Fourier 
modes  associated  with  the  velocity  field  Ui(x,y,z)  will  be  denoted  Ui(kx,ky,kz)  and  are 
defined  by 

Ui(kx,ky,kz)  = Ui(x,y,z)e  . (2-1) 

X 

Since  the  mixing  layer  is  homogeneous  in  two  directions  it  is  also  convenient  to  consider 
2D  Fourier  transforms.  Here  we  take  the  y direction  as  the  inhomogeneous  direction  and 
define  the  2D  Fourier  modes  Ui(kx,y,  kz)  by 

Ui{kx,y,kz)  = ^2  ui fol/, z)e“tkj-'Xj-,  (2.2) 

Xj. 

where  xj_  = (x,  z)  and  kx  = (kx,kz).  For  convenience,  we  will  also  adopt  the  following 
notations:  ux  = u,  uy  = v and  uz  = w.  When  the  flow  is  homogeneous  and  isotropic,  a 
common  way  to  initialize  the  modes  Ui(kx,  ky,kz)  is  to  fix  their  amplitudes  to  match  a 
given  energy  spectra  E(k)  and  assign  them  random  phases  in  such  a way  that  continuity 
is  enforced  (Rogallo  1981).  One  then  has, 

(\ui(kx,ky,kz)\2)  = A2(k),  with  E(k)  — 2nk2A2(k)  (2.3) 

and  k2  = k%  + k2  + k2.  For  the  case  at  hand,  we  can  adopt  a similar  strategy  but  consider 
instead  the  2D  spectra  in  each  plane  perpendicular  to  the  direction  of  inhomogeneity. 
Indeed,  in  those  planes  the  flow  is  assumed  to  be  homogeneous  and  isotropic.  We  thus 
initialize  our  flow  by  imposing  the  following  constraints  on  the  2D  Fourier  modes  (the 
assignment  of  the  random  phases  and  the  treatment  of  continuity  will  be  described  be- 
low), 

< \ui(kx,y,kz)\2  >=  B2(kx,y),  with  E(kx,y)  = nkxB2(kx,y)  (2.4) 

and  k\  = k2  + k2.  In  the  above  equation,  E(k± , y)  is  the  energy  spectra  of  velocity  field 
in  (x,  z)  plane.  The  arbitrary  part  remaining  is  the  choice  of  the  function  B2(kx,y).  For 
homogeneous  isotropic  flows,  it  is  trivial  to  relate  the  2D  amplitudes  to  the  3D  amplitudes 
(using  Parseval’s  theorem): 

B2{k±,y)  = J dkyA2(k2  + kx).  (2.5) 

Note  that  as  expected,  B2{kx,y)  is  independent  of  y for  homogeneous  flows.  In  the 
case  of  the  shear-free  mixing  layer  we  will  choose  an  amplitude  function  A(k)  for  each 
homogeneous  region  and  compute  the  corresponding  functions  B(kx,y)  using  (2.5).  If 
the  2D  amplitudes  functions  in  the  high  energy  and  low  energy  regions  are  respectively 
denoted  Bn(kx,y ) and  Bi(kx,y ),  we  then  define  the  complete  2D  amplitudes  function 
for  the  shear- free  mixing  layer  as, 

BML(k±,y)  = (1  - f(y))BL(kx,y)  + f(y)BH(kx,y).  (2.6) 

The  function  f(y)  is  equal  to  0 in  the  low  energy  region  and  equal  to  1 in  the  high 
energy  region;  inside  the  mixing  layers  it  varies  smoothly  from  0 to  1.  The  complete 
initialization  procedure  is  as  follows.  First  we  initialize  our  3D  Fourier  modes  using  the 
procedure  of  Rogallo  (1981).  The  3D  energy  spectra  used  here  is  taken  from  the  high- 
energy  homogeneous  region  of  the  mixing  layer.  This  ensures  that  the  Fourier  modes 
Ui(kx,  ky,  kz)  satisfy  the  continuity  equation.  The  3D  Fourier  modes  are  then  transformed 
to  Ui(kx,y,kz)  using  a ID  Fourier  transform.  At  this  point  their  2D  amplitudes  are 
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measured  and  the  modes  are  rescaled  in  order  to  match  the  prescribed  2D  amplitudes 
given  by  (2.6): 


Ui(kx,y,kz ) ->  u'i{kx,y , kz)  = Ui(kx,y,kz) * 


BML(k±,y) 
\ui (kx , y , A;j.)p 


The  are  then  transformed  back  into  3D  Fourier  modes,  u'i(kx,  ky,  kz).  By 

performing  the  transformation  (2.7),  one  of  course  destroys  the  continuity  property  of  the 
initial  field  and  it  has  to  be  recovered  by  projecting  the  ft((fcx,  ky,  kz)  onto  a divergence- 
free  field: 


kikj  \ . 


ky,  kz)  * (kx,  ky j kz ) — ^2  J ^ j{kxi  ky , (2.8) 

This  in  turn  partly  undoes  the  prescription  of  the  2D  amplitudes.  Fortunately,  by  iter- 
ating the  transformations  (2.7)  and  (2.8)  one  converges  to  a velocity  field  that  has  2D 
amplitudes  arbitrarily  close  to  BML(kj_,y)  and  which  satisfies  continuity.  At  this  stage, 
the  flow  still  has  random  phases.  In  order  to  correct  this  problem,  we  have  time  evolved 
the  flow  until  the  global  skewness  of  the  velocity  derivatives  reached  a converged  value. 
Between  each  time  step,  the  mixing  layer  was  rebuilt  using  (2.7)  and  (2.8)  in  order  to 
retain  the  desired  2D  amplitudes  profile.  After  that,  we  stopped  the  rebuilding  procedure 
and  let  the  flow  decay  freely. 


3.  DNS  results 

3.1.  Parameters  of  the  simulation 

The  choice  of  parameters  for  the  DNS  was  mainly  guided  by  the  following  considerations. 
In  order  to  have  an  experimental  reference  to  compare  with,  the  parameters  of  the  DNS 
have  been  chosen  to  match  as  closely  as  possible  those  from  the  Veeravalli  & Waxhaft 
(1989)  experiment  performed  using  the  3 : 1 perforated  plate.  Of  course,  since  numerical 
capabilities  are  not  unlimited,  some  compromises  had  to  be  made.  The  most  important 
restriction  in  the  present  study  is  the  ability  to  adequately  resolve  the  high-energy  region 
of  the  flow.  Given  this  constraint,  the  initial  3D  spectra  of  the  homogeneous  high-energy 
region  E{i(k)  was  chosen  to  match  the  spectra  measured  in  the  Comte-Bellot  & Corrsin 
(1971)  experiment  at  stage  1.  This  experimental  spectra  was  fitted  with  the  following 
function, 

EH{k)  = ?k—-e-bk*  (3.1) 

(k*  + q*)1+a  { ’ 

which  contains  several  parameters  a,  q,  6,  a and  (3.  This  fairly  complicated  function  has 
been  chosen  because  it  allows  an  easy  fit  of  various  properties  of  the  energy  spectrum.  For 
instance,  the  parameters  6 and  /?  can  be  used  for  characterizing  the  viscous  range  of  the 
spectrum.  The  parameters  q and  a determine  the  energy  peak  and  the  transition  between 
the  energy  containing  scales  and  the  viscous  range.  Finally,  a determines  the  total  energy. 
The  function  (3.1)  does  not  allow  one  to  derive  analytical  expressions  for  the  total  energy 
and  the  total  dissipation  in  terms  of  the  parameters  a,  q,  b,  a and  (3.  It  is  thus  not  possible 
to  express  these  parameters  in  terms  of  simple  global  experimental  data  and  we  have  not 
found  a systematic  procedure  for  prescribing  them.  It  has  however  been  observed  that 
the  following  set  of  parameters,  a = 10.6,  g = 1.5,  b = 0.02,  a = 1.233,/?  = 1.1,  allows  an 
almost  perfect  fit  to  the  Comte-Bellot  & Corrsin  spectrum.  Of  course,  the  value  of  some 
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Figure  2.  3D  contour  plot  of  the  energy  density  for  the  initial  velocity  field.  The  labels  indicated 
on  the  figure  correspond  to  grid-points.  The  data  field  was  sampled  every  four  grid-points  in 
order  to  produce  the  graph. 

of  these  parameters  might  depend  on  the  units  chosen  to  perform  the  simulation.  This 
is  not  really  an  issue  since  the  time  scale  and  the  length  scale  can  be  seen  as  entirely 
defined  by  the  computational  domain  size  lx  = 27T,  ly  = 47 r,  lz  = 2tt  and  by  the  viscosity, 
chosen  here  to  have  the  numerical  value  of  0.006.  In  the  Veeravalli  & Warhaft  (1989) 
experiment  the  ratio  of  energy  between  the  two  homogeneous  regions  is  6.27  while  the 
ratio  of  dissipation  is  7.28.  These  two  ratios  can  be  reproduced  well  by  choosing  the 
spectra  of  the  low  energy  homogeneous  region  to  be  of  the  same  form  as  (3.1)  but  with 
the  following  parameters:  a = 2.74,  q = 3.33,6  = 0.027,  a = 1.233,/?  = 1.1.  Furthermore, 
with  the  above  choices  of  parameters,  the  maxima  of  the  two  spectra  are  separated  by  a 
ratio  that  matches  the  inverse  ratio  of  the  initial  integral  length-scales  in  the  Veeravalli 
& Warhaft  (1989)  experiment  between  the  low-  and  high-energy  regions.  Accordingly, 
the  differences  in  typical  sizes  of  the  large-scale  structures  on  both  sides  of  the  mixing 
layer  are  also  reproduced.  From  these  definitions  it  is  possible  to  compute  numerically 
the  two  functions  Bi(k±,y)  and  Bj{(k±,y)  needed  in  (2.6)  using  (2.5).  The  initialization 
procedure  is  then  fully  defined  if  the  smoothing  function  f(y)  is  prescribed.  Here,  the 
following  choice  has  been  adopted: 

'0  if  0 < y < 5J„/24, 

l/2(sin(127r(v~iy-))  + 1)  if  5/y/24  < y < 7i„/24, 
f(y)  = 1 " if  7Zy/24  < y < 17ly/24,  (3.2) 

l/2(sin(127 r(y~y3))  + 1)  if  17/„/24  < y < 19J„/24, 
k0  if  19Zy/24  <y  < ly. 

With  this  smoothing  function,  the  high-energy  region  and  the  combined  low-energy 
regions  have  the  same  length.  Both  are  five  times  larger  than  each  of  the  two  mixing 
layers.  As  an  illustration,  the  contour  plot  of  the  initial  energy  density  is  shown  in  Fig.  2. 
As  far  as  the  numerics  are  concerned,  our  DNS  was  performed  using  a pseudo-spectral 
dealiased  code.  The  grid  resolution  adopted  consists  of  512  x 1024  x 512  points.  The 
higher  resolution  dimension  being  the  direction  of  inhomogeneity,  taken  here  to  be  y. 

3.2.  Kinetic  energy  diagnostics 

One  of  the  major  motivations  of  this  study  is  to  investigate  the  effect  of  inhomogeneities 
on  the  properties  of  the  turbulence.  Since  the  direction  of  inhomogeneity  is  along  y,  it 
is  convenient  to  present  statistics  obtained  by  averaging  over  the  x and  z directions. 
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Figure  3.  (a)  Energy  and  (b)  Dissipation  rate  profiles  across  the  mixing  layer  calculated  from 
(3.3)  and  (3.4); : t*  = 0; : f*  = 0.56; : t*  = 1.51. 

For  instance,  the  kinetic  energy  and  dissipation  rate  profiles  are  calculated  from  the 
expressions: 

E(y)  = JRWP  = -V  £ 5M*)I2  • (3-3) 

Z Tl/'rTiy  Z 

X,Z 

e(y)  = 2 v Sij  (x)5y  (x) , (3.4) 

where  S y (x)  = (dtu3  +djUi)/‘2.  Here,  and  in  the  rest  of  this  paper,  the  overbar  ~ denotes 
averaging  over  the  planes  perpendicular  to  the  direction  of  inhomogeneity.  Figures  3(a)  - 
3(b)  represent  the  profiles  of  the  kinetic  energy  and  the  dissipation  rate  at  three  different 
times  in  the  simulation.  Time  has  been  normalized  using  the  initial  eddy  turnover  time, 
t*  = te0/ko  where  t is  the  dimensional  time,  k0  is  the  initial  average  turbulent  kinetic 
energy  and  eo  is  the  initial  average  dissipation  rate.  The  Figs.  3 demonstrate  that  as 
the  decay  proceeds,  the  mixing-layer  widens  but  the  homogeneous  regions  remain  largely 
discernible.  Fig.  4(a)  shows  the  temporal  decay  of  the  average  energy  in  the  high-energy 
and  low-energy  homogeneous  regions.  Assuming  an  asymptotic  power-law  decay  E(t)  ~ 
tn,  a decay  exponent  of  n = —1.3  is  found  in  the  high-energy  homogeneous  region  while 
in  the  low-energy  region  the  decay  exponent  is  n = -1.1  (the  global  energy  decay  has  a 
decay  exponent  of  n = —1.3).  These  decay  rates  are  compatible  with  the  DNS  of  Briggs 
et  al.  (1996)  for  which  a kA  low  wave-vector  energy  spectrum  was  adopted  as  in  the 
present  simulation.  In  Fig.  4(b)  the  2D  spectra  defined  in  (2.4)  are  presented  to  confirm 
that  the  flow  is  sufficiently  resolved.  As  was  observed  in  Briggs  et  al.  (1996),  the  energy 
of  the  high  wave-numbers  decays  faster  with  time  than  the  energy  of  the  low  wave- 
numbers.  It  is  noted  that  for  the  discretization  used  and  the  times  considered,  the  energy 
peaks  remain  at  a constant  wave- vector  in  the  homogeneous  regions.  The  strong  influence 
from  the  homogeneous  layers  seems  to  induce  a shift  of  this  energy  peak  in  the  mixing 
layer  towards  higher  wave-vectors,  though  this  effect  remains  moderate  for  the  times 
considered. 


3.3.  Intermittency 

The  skewness  and  kurtosis  of  a velocity  component  u,  are  respectively  defined  as, 


KUi  - 


(«?) 


2)2 


(3.5) 
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Figure  4.  (a)  Energy  decays  of  the  high  (- ) and  low  ( ) energy  homogeneous  regions; 

( ) represents  the  global  energy  decay.  All  curves  have  been  normalized  using  their  initial 

value  of  fco.  (b)  2D  spectra  E(k±,y)  (defined  by  (2.4))  computed  inside  the  high  (H)  and  low 

(L)  homogeneous  regions  and  the  mixing-layer  (M)  for  three  different  times; : f*  = 0; 

: t*  = 0.56; : t*  = 1.51. 

For  homogeneous  isotropic  turbulence,  measures  of  these  quantities  show  that  they  are 
very  close  to  those  calculated  for  a Gaussian  signal,  i.e. , 5 = 0 and  K = 3. 

The  skewness  profile  of  v is  shown  in  Fig.  5(a).  The  y direction  has  been  normalized 
by  the  half- width  Zj/ 2 of  the  mixing  layer  and  centered  around  the  inflection  point  of 
the  variance  of  v.  Other  skewness  and  kurtosis  diagnostics  are  normalized  in  a similar 
fashion  (see  Veeravalli  & Warhaft  (1989)  for  the  details  of  this  normalization  procedure). 
The  skewness  profile  of  v exhibits  a sharp  deviation  from  the  Gaussian  value  around  the 
location  of  the  mixing  layer.  As  described  in  Veeravalli  & Warhaft  (1989),  this  behaviour 
is  attributed  to  the  intermittent  penetration  into  the  low-energy  region  of  structures 
originating  from  the  high-energy  region  (a  similar  penetration  of  structures  from  the 
low-energy  region  into  the  high-energy  region  is  certainly  also  happening  but  is  however 
a lot  less  frequent).  Agreement  with  experimental  data  from  Veeravalli  & Warhaft  (1989) 
is  very  good  both  in  terms  of  the  location  of  the  peak  and  its  amplitude.  For  symmetry 
reasons,  it  is  expected  that  the  skewness  of  u and  w should  remain  close  to  zero.  Up 
to  statistical  deviations,  this  is  confirmed  in  our  simulation  (although  not  illustrated  in 
this  paper).  Kurtosis  profiles  of  the  velocity  components  are  shown  in  Fig.  5(b,c,d)  and 
display  deviations  from  the  Gaussian  value  of  3 again  around  the  location  of  the  mixing 
layer.  In  the  Kurtosis  of  the  u,  we  observe  an  unexpectedly  high  peak  in  the  profile 
(compared  to  the  experimental  data).  At  this  point  we  attribute  this  issue  to  the  fact 
that  the  profile  was  computed  from  a single  realization  of  the  flow  (although  averaging 
in  the  (x,  z)  planes  was  performed).  We  also  note,  however,  an  initial  small  peak  in  this 
Kurtosis  component  present  in  the  initial  condition  (in  contrast  to  the  Kurtosis  of  the 
other  two  velocity  components)  but  it  is  not  clear  whether  or  not  it  was  amplified  during 
the  evolution  of  the  flow  and  by  what  mechanism.  The  other  two  Kurtosis  profiles  agree 
very  well  with  the  experimental  data  again  both  in  terms  of  the  location  of  the  peaks 
and  their  amplitudes.  Both  the  numerical  simulation  and  the  experiments  indicate  that 
the  deviations  from  the  Gaussian  values  for  5 and  K occur  on  the  low  energy  side  of  the 
inflection  points.  This  supports  the  idea  that  these  deviations  result  from  the  more  likely 
penetration  of  intermittent  structures  from  the  high-energy  region  into  the  low-energy 
region.  Finally,  it  must  be  stressed  that  the  initialization  procedure  presented  in  section 
2,  though  quite  sophisticated,  does  not  allow  imposition  of  the  initial  skewness  or  the 
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FIGURE  5.  (a)  Skewness  of  v;  (b)  Kurtosis  of  u;  (c)  Kurtosis  of  v;  (d)  Kurtosis  of  w; : 

f*  = 0; : t*  = 0.56; : f*  = 1.51;  o:  experimental  data  from  Veeravalli  & Warhaft 

(1989).  At  t*  = 0 the  numerical  curves  are  close  to  their  Gaussian  value  since  the  penetration 
mechanism  described  in  the  text  has  not  occurred  yet. 


kurtosis  profiles  in  the  region  of  the  mixing  layer  (Gaussian  values  are  in  fact  observed 
everywhere  for  the  initial  profiles  as  shown  in  Fig.  5 except  for  a small  deviation  for  Ku 
as  mentioned  above).  The  fact  that  the  DNS  later  reproduces  the  experimental  profiles 
observed  in  the  mixing  layer  indicate  that  the  transport  mechanisms  are  successfully 
resolved. 


4.  LES  of  the  shear-free  mixing-layer 

4.1.  Notations  and  conventions  - LES  model 

Starting  from  the  Navier-Stokes  equations,  one  obtains  the  LES  equations  (4.1)  by  ap- 
plying a filter,  here  denoted  ~ (since  our  code  is  spectral,  we  will  only  consider  spectral 
cut-offs  for  the  LES  filter): 

dtUi  + dj(ujUi)  = —dip  4-  i/Afij  - djfij.  (4.1) 

The  unknown  subgrid-scale  stress  (SGS)  tensor  fj,-  = upuj  — UiUj  needs  to  be  modelled 
in  terms  of  iq  in  order  to  close  (4.1).  In  this  work,  we  will  use  for  a model  proposed 
in  Wong  & Lilly  (1994)  and  further  studied  in  Carati  et  al.  (1995)  and  in  Dantinne  et  al. 
(1998).  This  model,  which  can  be  considered  as  a variant  of  the  dynamic  Smagorinski 
model  (Smagorinsky  1963;  Germano  et  al.  1991;  Lilly  1992;  Germano  1992),  has  been 
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shown  to  perform  very  well  in  the  context  of  homogeneous  isotropic  turbulence  and  its 
predictions  are  extremely  close  to  the  dynamic  Smagorinski  model.  The  advantage  of 
this  model  rests  upon  the  ease  with  which  its  dynamic  version  can  be  implemented.  The 
definition  of  the  model  is  the  following: 

Dj  -^fkkSij  = -2CA*Sij,  (4.2) 

where  — \{dtUj  + djUt)  is  the  resolved  strain  tensor  and  A is  the  LES  filter  width. 
The  dimensional  parameter  C is  evaluated  by  introducing  a second  (coarser)  filter  — 
(the  test  filter)  and  using  the  dynamic  procedure: 

C = — r ^ X < A'--  (4.3) 

2(A4/3  - A4/3)  < SijSij  > 

where  Ly  = Uiiij  — utUj  is  the  Leonard  tensor  (note  that  we  have  systematically  used  the 
property  ~ = rr.  valid  for  spectral  cut-offs).  As  in  the  dynamic  Smagorinski  model,  the 
only  free  parameter  available  is  the  ratio  of  the  spectral  cut-offs:  A/A.  In  the  following 
discussion  this  ratio  will  be  assumed  to  be  equal  to  2.  For  homogeneous  isotropic  turbu- 
lence, the  averages  (■  • • ) in  (4.3)  are  obtained  by  averaging  over  the  whole  computational 
domain.  The  idea  is  of  course  that,  since  turbulence  is  homogeneous,  the  constant  C 
should  be  statistically  independent  of  the  position.  For  inhomogeneous  flows  in  one  di- 
rection, like  the  shear-free  mixing  layer  or  the  channel  flow,  dependence  on  the  direction 
of  inhomogeneity  is  introduced  by  averaging  quantities  only  over  the  other  two  homoge- 
neous directions.  This  is  justified  only  if  the  flow  is  not  too  inhomogeneous.  The  dynamic 
coefficient  C then  depends  explicitly  on  the  inhomogeneous  direction:  C — C(y).  Our 
LES  condition  was  obtained  by  filtering  (with  a sharp  spectral  cut-off)  the  initial  DNS 
field  down  to  32  x 64  x 32  modes.  Thus,  only  of  the  gridpoints  are  retained  in  each 
direction  for  the  LES,  meaning  that  there  is  one  LES  grid  point  for  about  4 000  DNS 
grid  points.  The  box-size  is  unchanged  and  remains  2tt  x x 27t.  In  our  study  we  have 
also  included  a run  obtained  at  LES  resolution  but  with  no  subgrid-scale  stress  model  to 
emphasize  the  effect  of  the  model  in  the  LES  simulation.  For  comparison  we  have  filtered 
the  DNS  fields  stored  during  the  simulation  down  to  32  x 64  x 32  modes. 

4.2.  Comparison  of  the  filtered  DNS  and  the  LES 

Figure  6(a)  represents  the  temporal  evolution  of  the  normalized  global  kinetic  energy 
E/Eq.  From  the  graph  it  is  clear  that  the  simulation  with  LES  model  does  a very  good 
job  in  reproducing  this  diagnostic,  while  the  run  with  no  model  does  not.  Likewise,  Figs. 
6(b)  and  6(c)  show  that  the  2D  spectra  with  LES  modelling  axe  in  good  agreement  with 
DNS  data.  2D  spectra  gathered  from  the  run  with  no  model  exhibit  a strong  pile  up  of 
energy  in  the  high  wave-number  side  of  spectra,  indicating  that  the  flow  is  not  adequately 
resolved  in  that  case. 

To  further  illustrate  the  decay  of  the  kinetic  energy,  we  display  in  Fig.  7(a)  and  7(b) 
the  profiles  of  the  kinetic  energy  at  two  different  times.  This  is  of  course  a more  sensitive 
diagnostic  since  it  retains  information  about  the  inhomogeneity  of  the  flow.  Both  figures 
again  reveal  a good  agreement  between  the  DNS  and  LES  runs  and  a poor  behaviour 
of  the  run  without  modelling.  Variance  profiles  of  u,  v and  w have  been  examined  (not 
shown)  and  again  the  LES  matches  the  DNS  very  nicely,  whereas  the  no  model  simulation 
performs  poorly.  For  our  assessment  of  LES  we  have  retained  the  same  intermittency 
diagnostics  described  earlier:  the  skewness  and  the  kurtosis  of  the  velocity  components. 


E/Eo 


Figure  6.  (a)  Evolution  with  time  of  the  normalized  global  kinetic  energy  E/Eo\  (b)  and  (c) 
2D  spectra  E(kj_,y)  (defined  by  (2.4))  computed  inside  the  high  (H)  and  low  (L)  homogeneous 
regions  for  time  t"  = 1.51.  In  both  figures, : filtered  DNS; : LES; : no  model. 


Figure  7.  Energy  profile  across  the  mixing-layer  calculated  from  (3.3)  at  times  (a)  t * = 0.56, 
(b)  t*  = 1.51.  In  both  figures, : filtered  DNS; : LES; : no  model. 


A sample  of  those  quantities  is  displayed  in  Fig.  8.  Once  again,  the  LES  run  produces 
results  which  compare  well  with  the  filtered  DNS  data.  Surprising  at  first,  the  run  without 
subgrid-scale  stress  model  also  performs  quite  well.  However,  since  the  intermittency  is 
attributed  to  large-scale  structures  penetrating  the  low  energy  region  from  the  high 
energy  region  this  is  to  be  expected.  Indeed,  looking  at  the  spectra  displayed  in  Fig.  6(c), 
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Figure  8.  (a)  Skewness  of  v\  (b)  Kurtosis  of  u;  (c)  Kurtosis  of  v;  (d)  Kurtosis  of  w.  All  curves 
computed  at  t*  = 1.51. : filtered  DNS; : LES; : no  model. 


we  see  that  the  run  without  subgrid-scale  stress  model  is  not  ill-behaved  for  modes  at 
the  lowest  wave-numbers. 


5.  Conclusions 

This  study  of  the  shear-free  mixing  layer  reveals  that  DNS  is  capable  of  reproducing 
most  of  the  aspects  of  the  experimental  database  that  have  been  produced  by  Veeravalli 
& Warhaft  (1989).  Energy,  dissipation  rate  and  velocity  variance  profiles  are  accurately 
reproduced.  Also,  departure  from  Gaussianity  inside  the  mixing  layer  revealed  by  the 
measurement  of  the  skewness  and  the  kurtosis  of  the  velocity  field  are  also  in  excellent 
agreement  with  the  experimental  observation.  This  is  particularly  satisfactory  since  inside 
the  mixing  layer  the  initial  values  of  these  quantities  cannot  be  prescribed  by  the  initial- 
ization procedure  and  are  entirely  produced  at  later  times  by  the  transfer  mechanisms 
simulated  by  the  DNS.  It  is  also  observed  that,  in  the  absence  of  shear,  no  significant 
anisotropy  is  observed  in  this  flow.  This  is  consistent  with  previous  numerical  results  and 
strongly  supports  the  use  of  eddy-viscosity  type  models  for  the  LES  of  such  flows.  As  a 
consequence,  it  is  not  surprising  that  comparisons  between  the  predictions  of  LES  using 
such  eddy-viscosity  models  and  DNS  show  very  good  agreement. 
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